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Chapter  1 

SUMMARY  OF  THE  PROJECT 


1.1  RESEARCH  OBJECTIVES 

The  general  objective  of  this  research  is  to  study  the  mechanics  of 
granular  material  from  a  micro-structural  point  of  view.  The  specific 
objectives  of  this  research  are  as  follows: 

1.  development  of  a  mathematical  model  to  obtain  the  stress-strain 
relationship  for  a  granular  material  explicitly  considering  the  micro¬ 
structure  characteristics  of  the  packing, 

2.  study  of  the  force-displacement  law  at  a  contact  between  two  non- 
conforming  bodies  and  its  implementation  in  the  stress-strain  relationship, 

3.  development  of  the  mathematical  characterization  of  the 
micro-structure  of  a  granular  assembly  and  investigation  of  its  effect  on 
the  stiffness  characteristics  of  the  packings, 

4.  experimental  study  on  rod  assemblies  in  a  directional  shear  box  to 
investigate  the  effects  of  packing  structure  and  loading  conditions  on  the 
mechanical  behavior  of  granular  assemblies,  and 

5.  verification  of  the  theoretical  model  by  comparing  the  theoretical 
results  with  those  from  tests  on  rod  assemblies  in  directional  shear  box  and 
other  available  tests  in  the  literature. 


1 . 2  ACCOMPLISHMENTS 
1.2.1  Theory 

A  general  mathematical  model  is  developed  to  evaluate  the  stress-strain 
behavior  of  granular  packings  based  on  micro-structural  considerations.  Two 


approaches  are  followed,  namely,  (a)  mobilized  plane  approach  and  (b) 
particulate  approach. 

(a)  The  mobilized  plane  model  assumes  that  a  large  number  of  mobilized 
planes  are  developed  in  a  soil  element  when  subject  to  loading.  The 
overall  mechanical  behavior  of  a  packing  can  be  obtained  based  on  the 
movement  of  the  mobilized  planes.  Work  on  this  approach  is  given  in 
publications  1,  2,  4  and  7  shown  in  section  1.3.3. 

(b)  The  particulate  model  explicitly  incorporates  the  geometrical 
characteristics  of  the  packing,  such  as,  the  void  ratio,  the 
coordination  number,  and  the  contact  normal  distribution.  The  work 
accomplished  in  this  approach  is  as  follows: 

(b.l)  The  formulation  of  the  stress-strain  relationship  for  regular 

packings  of  equal  size  particles  is  given  in  publications  5  and 
12  in  section  1.3.3. 

(b.2)  The  formulation  of  the  stress-strain  relationship  for  regular 
packings  of  multi  size  particles  is  given  in  publication  8  in 
section  1.3.3. 

(b.3)  The  formulation  of  the  stress-strain  relationship  for  random 

packings  of  multi  size  particles  is  given  in  publications  3  and 
10  in  section  1.3.3. 

The  theoretical  development  for  the  particulate  model  is  discussed  in 
Chapter  2  of  this  report. 


1.2.2  Contact  Law 

The  contact  force-displacement  lav  for  two  non-conforming  bodies  in 
contact  is  incorporated  in  the  stress-strain  formulation  for  the  particulate 
approach.  The  work  accomplished  is  as  follows: 
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(a)  Closed-form  solutions  of  the  constitutive  matrix  for  linear 
contact  lav  are  given  in  publications  6,  9  and  10  in  section  1.3.3. 

(b)  Solutions  for  non-linear  Her tzian-Mindlin  contact  lav  are  given  in 
publications  6,  9,  10  and  11  in  section  1.3.3. 

The  contact  lav  is  discussed  in  Chapter  3  and  the  solutions  are  shovn 
in  Chapter  4  and  6  of  this  report. 


1.2.3  Fabric  and  its  Effect  on  the  Mechanical  Properties 

Mathematical  characterization  of  the  micro-structure  of  a  granular 
assembly  is  developed.  This  is  incorporated  in  the  stress-strain 
formulation  for  the  particulate  approach  to  study  the  effect  of  packing 
structure  on  the  stiffness  properties.  The  vork  accomplished  is  as  follovs: 

(a)  The  formulation  for  representing  contact  normal  distribution  and 
its  relationship  vith  the  fabric  tensor  is  given  in  publications  9  and 
10  in  section  1.3.3. 

(b)  The  stiffness  properties  of  packings  vith  anisotropic  micro¬ 
structure  are  given  in  publications  9,  10  and  11  in  section  1.3.3. 

The  discussion  on  fabric  and  its  effect  on  stiffness  properties  of  a 

packing  is  given  in  Chapter  4. 


1.2.4  Experiments  on  Rod  Assemblies 

Experiments  are  conducted  on  assemblies  of  cylindrical  rods  to 
investigate  the  deformation  behavior.  The  vork  accomplished  is  as  follovs: 

(a)  Results  on  directional  shear  box  are  given  in  publications  5,  8 
and  13  in  section  1.3.3. 

(b)  Results  on  model  footing  tests  are  given  in  publication  14  in 
section  1.3.3. 
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The  experimental  results  on  the  directional  shear  cell  is  presented  in 
Chapter  5  of  this  report. 

1.2.5  Verification  of  the  Theoretical  Model 

The  theoretical  model  is  verified  by  comparing  the  theoretical  results 
with  those  from  the  tests  in  directional  shear  box  and  other  available  test 
results  in  the  literature.  The  work  accomplished  is  as  follows: 

(a)  Comparison  with  results  from  tests  in  directional  shear  box  are 
given  in  publications  5  and  8  in  section  1.3.3. 

(b)  Comparisons  with  results  on  assemblies  of  glass  balls  in  resonant 
column  device  and  sands  under  low  amplitude  waves  are  given  in 
publication  11  in  section  1.3.3. 

The  comparisons  are  discussed  in  Chapter  6  of  this  report. 
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-  "Behavior  of  Plastic  Flow  for  Sand," 

2.  Engineering  Mechanics  Division,  Specialty  Conference,  ASCE,  Buffalo, 
New  York,  May  1987. 

-  "A  Micro-structural  Approach  to  Constitutive  Modelling  for 
Granular  Soils," 

3.  International  Workshop  on  Constitutive  Equations  for  Granular  Non- 
cohesive  Soils,  Cleveland,  Ohio,  July  22-24,  1987. 

-  "Micromechanical  Approach  to  Constitutive  Modelling  for  Granular 
Sand  Based  on  Slip  Mechanism," 

4.  U.S.  Japan  Seminar  on  the  Micromechanics  of  Granular  Materials, 
Sendai,  Japan,  Oct.  26  -  30,  1987. 

-  "Micromechanical  Modelling  of  Constitutive  Equation  for  Granular 
Material, " 

5.  Engineering  Mechanics  Division,  Specialty  Conference,  ASCE, 
Blacksburg,  Virginia,  May  1988. 

-  "Deformation  Analysis  for  Granular  Material," 

6.  Engineering  Mechanics  Division,  Specialty  Conference,  ASCE, 
Blacksburg,  Virginia,  May  1988. 

-  "Cyclic  Stress-strain  Modelling  for  Granular  Soils  with  Micro- 
structural  Consideration," 
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Chapter  2 

THEORETICAL  DEVELOPMENT 


The  stress-strain  modelling  of  discrete  systems,  such  as  granular 
packings,  can  only  be  acheived  by  considering  the  continuum  field  quantities 
as  averages  of  their  discrete  analogues  (Chang  1987,  Chang  and  Misra  1988b). 
The  stress,  a  continuum  field  variable,  is  related  to  the  discrete  contact 
forces  in  a  granular  media.  Similarly,  the  strain  in  a  granular  media  is 
related  to  the  contact  displacements. 


$ 


I 

g 

9 


v 

ft 


V 

,  v 

n\ 


I* 


JO 
!  * 


% 

V 

» 


2.1  STRESSES  AND  CONTACT  FORCES 

Assembly  of  spherical  particles  can  be  divided  into  polyhedral  sub¬ 
volumes  defined  as  'Voronoi  cells'  such  that  the  space  occupied  by  the 
assembly  is  completely  filled  (see  Fig.  2.1).  This  polyhedral  cell  is 
chosen  to  contain  a  particle  and  the  void  space  closest  to  that  particle 


(Finney  1970).  The  average  stress  tensor  within  the  nth  (see  Fig.  2.1) 


(2.1) 


'Voronoi  cell'  is  given  by 

a? .  =  — f  a.  .  dv 
AVn  J  1J 

where  AVn  is  the  volume  of  the  nth  'Voronoi  cell'.  From  Eq.  2.1  and  using 
the  theorem  of  stress  means  (Truesdell  and  Toupin  1960),  the  average  stress 
tensor  for  the  'Voronoi  cell'  is  expressed  in  terms  of  surface  tractions  on 
the  'Voronoi  cell'  as  follows 

a?.  =  —  f  t.  r.  ds  (2.2) 

U  AVn  J  1  3 


where  t^  is  the  traction  vector  on  the  surface  of  the  nth  sub-volume  at  the 


location  given  by  the  position  vector  r  ^  with  the  origin  chosen  at  the 


centre  of  the  sphere  contained  in  the  nth  sub-volume.  In  a  packing  of 
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spheres,  instead  of  continuous  surface  tractions  there  are  discrete  contact 

points  on  the  surface  of  the  nth  sub-volume  where  the  forces  act  (see  Fig. 

2.1).  Therefore,  for  this  system,  Eq.  2.2  reads  as 

n.  =  —  T  tnm  rnm  (2.3) 

™  1  J 


-n 

<3. 

1J 


AV  m 
nm 


where  ti  is  the  contact  force  at  the  mth  contact  of  the  nth  particle  and 
r?m  represents  the  position  vector  of  that  contact. 

For  continuum  modelling  of  packings  of  spheres,  it  is  necessary  to 
define  the  stress  for  a  representative  volume  which  best  represents  the 
packing.  This  representative  volume  comprises  ot  many  'Voronoi  cells'  of 
the  packing.  Clearly,  the  stress  field  in  an  assembly  of  spheres,  due  to 
its  discrete  nature,  is  heterogeneous,  that  is,  the  average  stress  of 
the  'Voronoi  cell'  varies  from  one  cell  to  another.  Therefore,  in  this 
case,  a  'mean'  stress  tensor  for  a  chosen  representative  volume  may  be 
defined.  The  'mean'  stress  tensor,  for  this  representative  volume,  is 
defined  as  an  average  of  the  stress  tensors  of  its  sub-volumes  ('Voronoi 
cell')  as  follows, 

<°ij>  E  AV°  (2-M 

where  is  the  average  stress  tensor  for  the  sub-volumes  AVn  and  the 


volume  V  is  given  by 

v  =  jr  avn 


(2.5) 


Further,  substituting  Eq.  2.3  into  Eq.  2.4,  the  'mean'  stress  tensor  of  the 
medium  can  be  written  as 

.  .  1  t—  f  .nm  nm  . 

<”(ij)>  -  »  E  E  <(1  -j,  <2-6) 

J  n  m  x  J ' 

where  ()  represents  the  symmetric  part.  The  skew  symmetric  part,  given  by, 


<a 


.  1  t—  t—  nm  nm 

ill5  ■  v  E  E  t(i  r., 


(2.7) 


represents  the  'mean'  torque  of  the  representative  volume  contributed  by  the 
torque  for  each  sub-volume  caused  by  the  contact  force  t^m  acting  at  r^m. 

In  absence  of  any  applied  torque  or  couple  stresses,  the  skew  symmetric  part 
is  expected  to  be  zero  in  order  to  satisfy  Cauchy's  law  of  balance  of 
angular  momentum  (Truesdell  and  Toupin  1960). 

It  is  noted  that  representation  of  stress  tensor  in  granular  media, 
similar  to  Eq.  2.6,  have  also  been  developed  by  Chrisstof fersen  et.  al. 

1981,  Rothenburg  and  Selvadurai  1981,  and  Kishino  1978  using  different 
arguments. 

2.2  KINEMATICS  OF  GRANULAR  MEDIA 

2.2.1  Relative  displacement  betveen  two  contacting  spheres 

Under  an  increment  of  loading,  the  spheres  comprising  the  granular 
assembly  displace  from  their  original  position  leading  to  the  deformation  of 
the  assembly.  Each  particle  in  the  assembly,  in  general,  has  two  modes  of 
movement,  namely,  translational  movement  denoted  by  vector  u  and  rotational 
movement  denoted  by  o>.  These  movements  are  accompanied  by  the  interactions 
at  the  contact  between  two  particles.  Assuming  the  particles  to  be  rigid, 
the  relative  displacement  at  a  contact  betveen  two  particles,  n  and  p  (see 
Fig.  2.2),  can  be  expressed  by 

5  =  (u11-  up)  +  («n  x  r11)  +  (wP  x  rp)  (2.8) 

where  the  first  term  on  the  right  hand  side  of  Eq.  2.8  is  contributed  by  the 
translation,  while  the  other  terms  are  from  particle  spin.  Further,  the 
relative  displacement  S  can  be  expressed  as 

6  =  6nm  +  Spm  (2.9) 

where  6nm  and  are  the  relat  ive  displacements  associated  with  particles  n 
and  p  respectively,  and  superscript  m  represents  the  contact  between  n  and 
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p.  The  relative  displacement  6  has,  in  general,  tvo  components,  namely,  the 
normal  relative  displacement  and  the  shear  relative  displacement. 

2.2.2  Strain  tensor 

Similar  to  the  stress  field,  the  strain  field  in  a  packing  of  spheres 
is,  in  general,  inhomogeneous  and  discontinuous.  Analogous  to  the  ’mean’ 
stress  tensor,  a  'mean'  strain  tensor  can  be  established  for  the  chosen 
representative  volume.  The  'mean'  strain  tensor  for  the  representative 
volume  is  defined  from  energy  equivalence,  which  states,  in  general  terms, 
that  the  work  done  expressed  in  terms  of  the  'mean'  stress  and  'mean'  strain 
tensor  is  same  as  the  work  done  expressed  in  terms  of  the  contact  forces  and 
the  relative  displacements  at  the  contacts. 

The  work  done,  for  the  representative  volume  V,  expressed  in  terms  of 
the  'mean'  stress  and  the  'mean'  strain  is  given  by 

dW  =  V  <ff(ij)>  <feij>  (2.10) 

where  <e^>  Is  the  increment  of  the  symmetric  'mean'  strain  tensor.  By 
definition,  the  increment  of  symmetric  'mean'  strain  tensor  is  related  to 
the  increment  of  displacement  gradient  <D^>  and  the  increment  of  the 
rotational  tensor  <2^>  as  follows, 

<e. .>  =  <d. .>  -  <a. .> 

ij  ij  ij 


(2.11) 

Needless  to  say  that,  as  the  stress  tensor  is  symmetric,  the  rotational 


tensor  <Q^>  (i.e.,  the  anti-symmetric  part  of  the  displacement  gradient) 
does  not  contribute  to  the  work  done,  dW. 

The  work  done  expressed  in  terms  of  the  contact  forces  and  the 
incremental  relative  displacements  at  the  contacts  for  the  same 
representative  volume,  is  given  by 


1  .* 
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(2.12) 


dW  =  £  I!  t?m 

n  n  1  1 

Equating  Eq.  2.10  and  2.12  and  substituting  the  expression  for 
<°(ij)>  from  Eq.  2.6,  one  obtains  the  relationship  between  the  relative 
displacement  and  the  'mean'  strain  tensor 

(2.13) 

Using  Eq.  2.11,  Eq.  2.13  can  be  further  written  in  terms  of  the  displacement 
gradient  tensor  and  the  rotation  tensor  as  follows 

(2.14) 


6™ 

l 


nm  . 

r .  <e.  .> 

J  iJ 


6nm  =  rnm  <D. .> 

i 


r?m  <2. .> 

J  ij 


J  iJ 

2.2.3  Particle  spin  and  rotation  tensor 


In  Eq.  2.14,  6?m  represents  the  relative  displacement  at  the  mth 
contact  of  the  nth  ball  with  respect  to  its  contacting  neighbor  p 
considering  that  the  pth  ball  does  not  move.  Further,  by  definition,  the 
relative  displacement  of  the  nth  and  the  pth  balls,  compatible  with  the 


displacement  gradient  field  D^,  is 

n  p  nm 

u  -  =  r 


<D, ,  >  -  r?m  <D. .> 


(2.15) 


j  -if  -j  -i  y 
Thus  comparing  Eqs.  2.8,  2.14  and  2.15,  the  following  relationship  is 

obtained  between  the  particle  spin  and  the  rotation  tensor, 

r11  x  wn  +  rp  x  up  =  [  r?m  +  r^m  j  <®ij>  (2.16) 

Eq.  2.16  represents  an  underlying  assumption  that  the  mean  particle  spin  is 

same  as  the  'mean'  rigid  body  rotation  2^  of  the  representative  volume. 

2.3  RELATION  BETWEEN  INCREMENTAL  STRESS  AND  STRAIN 
2.3.1  Local  constitutive  lav 

A  local  constitutive  law  relating  the  incremental  contact  force  and 
the  relative  displacement  at  the  contact  is  defined.  For  convenience,  the 
local  constitutive  law  is  defined  in  a  local  cartesian  coordinate  system 


formed  by  the  base  vectors  n,  s  and  t.  Where  n  is  the  unit  normal  vector  at 
the  mth  contact  and  s  and  t  are  chosen  arbitrarily  in  genera  1  (Fig.  2.3). 

Let  T?™  denote  the  transformation  between  the  local  coordinate  system  and 
the  global  coordinate  system  such  that 

d.  =T.k&k,  and  f.  =  T.k  tk  (2.17) 

where  f?m  and  are  the  force  increments  at  the  mth  contact  of  the  nth 

particle  in  the  local  and  global  coordinate  system  respectively,  dnm  and  6^m 

J  it 

are  the  relative  displacements  increment  at  the  mth  contact  of  the  nth 
particle  in  the  local  and  global  coordinate  system  respectively.  The 
transformation  tensor  T-k,  from  the  global  to  the  local  coordinate  system, 
is  given  by 

T  _  «(j  )  /  T  1  O  \ 


T..  =*  e. 
jk  k 


(2.18) 


where  eVJ  =  {n,  s,  t)  is  a  set  of  basis  vector  forming  the  local  coordinate 
system,  such  that  e^*n,  e^=s,  and  e^  =  t.  Vector  n  is  defined  as 


n  =  (  cosy,  siny  cosg,  siny  sin3) 


(2.19) 


s  is  given  by 


s  =  (  -siny,  cosy  cos0,  cosy  sin(3) 


(2.20) 


and  t  by 


t  =  (  0,  -sing,  cos3) 


(2.21) 


where  y  and  3  are  defined  in  Fig.  2.3. 

If  the  local  constitutive  law  is  defined  in  the  local  coordinate  system  as 


follows, 


rnm  . , n rn  .nin 

t .  =  K. .  .  a  . 

i  ij  j 

.  .1 _ nm  . 


(2.22) 


where  is  the  stiffness  tensor.  Then,  in  the  global  coordinate  system, 


the  local  constitutive  law  reads  as  follows, 
.nm  _nm  „nm  _nm  .nm 

1  1  ~  1  ^  •  lx,,  i  .  ■  O. 

1  li  ij  jk  k 

2.3.2  Constitutive  law  for  the  packing 


(2.23) 


Fig-  2.3  l.ocal  coord 


Using  the  local  constitutive  law  defined  above,  the  overall 
constitutive  relationship  of  a  packing  of  spheres  can  be  derived. 


Substituting  for  6^  from  Eq.  2.13  and  multiplying  by  the  position  vector 


r  on  both  sides,  Eq.  2.23  becomes 

q 

nm  .nm  nm  Tnm  ..nm  Tnm  nm  . . 
r  t  >  s  r  i  i  >  ,  I,!  r  >  c,  / 


(2.24) 


A.  *.  i  ,  ,  IN  .  .  a.  ,  .  i.  NO,/'  V  ^  7 

q  1  q  li  ij  jk  p  kp  v 

Summing  up  on  both  sides  of  Eq.  2.24  and  dividing  by  volume  V,  the 

relationship  between  incremental  stress  and  the  incremental  strain  is 
obtained  as  follows 


<a.  >  =  C,  .  <e.  > 

lq  lqkp  kp 

where 


1  r—  r~  nm  _nm  „nm  Tnm  nm 
V  *—  fr  rq  li  ij  jk  rp 


(2.25) 


(2.26) 


and  the  incremental  stress  is  defined  as 


<V?  E  E‘I 


nm  nm 

1  rq 


(2.27) 


It  can  be  shown  that  the  incremental  stress  and  the  incremental  strain  in 
Eq.  2.26  are  frame  indifferent  quantities,  thus  satisfying  the  principle  of 
material  objectivity.  For  example,  let  Q  be  any  orthogonal  transformation 
tensor  from  coordinate  system  X  — »  X  ,  such  that 


.nm  -  ,nm  nm  nm  „ 
t  =  Q  .  t .  ;  r.  =r.  Q ., 
a  ai  l  b  j  jb 


(2.28) 


where  fnm  and  rPm  are  in  the  transformed  coordinate  system.  Then  the 
a  b 

transformed  incremental  stress  is  given  by 


.  *  t  1  r—  r— 

<  b>  =  Q  <  b>  Q  =  ^  H  t 


•k  k 

nm  nm 
a  rb 


(2.29) 


implying  that  the  incremental  stress  defined  in  Eq.  2.27  is  a  frame 
indifferent  quantity  under  any  orthogonal  transformation.  Similarly  the 


incremental  strain  can  also  be  shown  to  be  a  frame  independent  quantity. 


2.4  INCREMENT  OF  'MEAN'  STRESS  TENSOR  AND  INCREMENTAL  STRESS 


The  increment  of  the  'mean'  stress  tensor  is  defined  as 

<Ac.  .>  =  <a.  .>  -  <(j?.>  (2.30) 

ij  13  ij  v  ' 

where  <<r„>  and  are  the  'mean'  stress  tensors  in  the  deformed  and  the 

original  configurations  respectively.  Clearly,  the  incremental  stress 

<a..>,  defined  in  Eq.  2.27,  is  not  same  as  the  increment  of  the  'mean' 
ij 

stress  tensor.  In  order  to  express  the  constitutive  equation  in  terms  of 
the  increment  of  'mean'  stress  tensor,  the  relation  between  the  incremental 
stress  and  the  increment  of  the  'mean'  stress  tensors  needs  to  be 
established. 

From  Eq.  2.6,  the  'mean'  stress  tensors  in  the  deformed  and  the 
original  configurations  respectively  are  given  by 

•  J  L  Eti  r  and  <»°  >  .  -L  £  £  t°  (2.31) 


<’u> 


n  m  J  V  n  m 

The  superscript  nm  has  been  dropped  in  Eq.  2.31  for  clarity.  The  contact 

force  tj  and  the  position  vector  r^  in  the  deformed  configuration  can  be 

expressed  in  terms  of  the  contact  force  t?  and  the  position  vector  r°  in  the 

original  configuration  as  follows 

t.  =  t?  +  t. 

111 


(2.32) 

(2.33) 

where  the  dot  refers  to  the  incremental  change.  The  incremental  change  in 


r  ,  =  r  .  +  r  . 
J  J  J 


the  position  vector  r j ,  by  definition,  is  given  by 

r.  =  <D..>  r.  (2.34) 

where  <D^>  is  the  Eulerian  displacement  gradient.  Further,  the  Tr(D)  is 
defined  as 


Tr(D)  =  <Dnn>  *  1  -  ~ 


(2.35) 


18 


Substituting  Eqs.  2.31,  2.32,  and  2.33  into  Eq .  2.30,  using  Eqs.  2.34 

and  2.35  and  neglecting  all  the  quadratic  terms,  the  relation  between  the 

increment  of  the  'mean'  stress  tensor  <A<r.  .>  and  the  incremental  stress 

ij 

tensor  <c..>  is  obtained  as  follows 
ij 


<&cr.  .>  =  <cr.  . >  +  <0.  >  <D  .>  -  <0.  .>  <D  > 
lj  lj  lm  mj  ij  nn 


(2.36) 

Eq.  2.36  represents  the  familiar  relationship  between  the  increment  of 


Cauchy  stress  <Ao. .>  and  the  increment  of  first  Piola-Kirchof f  stress  <o.  .> 
ij  lj 

(Truesdell  and  Noll  1965).  If  the  symmetry  of  the  incremental  stress  tensor 

<a.  .>  is  rendered  explicit,  Eq.  2.36  becomes 
ij 

<d.  .>  =  <6cr.  . >  -  \  [  <0.  >  <D  .>  +■  <0.  >  <D  .>  1+  0.  .  <D  >  (2.37) 

ij  ij  2  L  im  mj  jm  mi  J  ij  nn 

Thus  the  constitutive  equation  can  be  written  in  terms  of  the  increment  of 
the  Cauchy  stress  tensor  as  follows 

<&<y.  .>  =  C....  <£,,>+  <a.  >  <D  .>  +  <0.  >  <D  .>  1-  0.  .  <D  >  (2.38) 

lj  ljkl  Tel  2  L  im  mj  jm  mi  J  ij  nn 

For  the  case  of  small  deformations  the  increment  of  the  Cauchy  stress  can  be 
assumed  to  be  same  as  increment  of  the  first  Piola-Kirchof f  stress,  that  is 
the  second  and  the  third  terms  on  the  right  side  of  Eq.  2.37  are  negligible 
compared  to  the  first  term.  In  that  case,  Eq.  ^38  reads 

<AV  =  Ci jkl  <£kl>  (2-39> 

where  same  as  defined  in  Eq.  2.26. 

2.5  CONSTITUTIVE  EQUATION  FOR  A  CONTINUOUS  SYSTEM 

The  constitutive  equation  discussed  thus  far  is  defined  for  the  case 
of  a  discrete  distribution  of  the  contact  position  vectors.  However,  for  a 
random  packing,  as  the  number  of  contacts  becomes  very  large,  the  contact 
position  vector  distribution  can  be  reasonably  treated  as  continuous.  A 
distribution  function  l( Q)  of  the  contact  position  vectors  can  be  defined 


such  that 


Thus  the  summation  in  Eq.  2.26  can  be  replaced  by  integration  and  the 
constitutive  tensor  becomes 

Clqpk  ■  f  Is  Tli<2>  V0’  Tjk(0)  rp<!!)  «s> da  <2-41) 

where  N  is  the  total  number  of  contacts  in  volume  V,  dSJ  =  siny  dy  d0  is  the 
differential  solid  angle  associated  with  a  unit  sphere  with  0<y<it  and  0<3<2 
and  £(S2)  represents  the  distribution  of  the  contact  position  vectors  in 
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Chapter  3 

CONTACT  LAW  UNDER  CYCLIC  LOADING 

THe  force-displacement  behavior  at  a  contact  significantly  effects 
the  stress-strain  behavior  of  a  granular  packing.  The  force-displacement 
law  of  two  non-conforming  bodies  in  contact  can  be  obtained  from  the  Hertz- 
Mindlin  theory  of  frictional  contacts.  The  normal  and  the  tangential 
stiffnesses  based  on  this  theory  are  presented.  The  assumptions  of  this 
theory  and  their  consequences  on  the  contact  law  are  discussed. 

3.1  NORMAL  CONTACT  OF  TWO  NON-CONFORMING  SOLIDS 

3.1.1  Hertz  theory 

For  two  spheres  of  radii  R^  and  elastic  properties  Gj  and  v,. 
(i=l,2)  the  classical  theory  of  Hertz  relates  the  normal  force  N  to  the 
relative  approach  a  of  the  centres  of  the  two  contiguous  spheres  and  the 
radius  a  of  the  circular  contact  area  by  the  following  equations  (Johnson 


(3.1) 

(3.2) 

(3.3) 

(3. A) 
(3.5) 


The  incremental  stiffness  is  defined  as 


AN 


=  2E'a 


(3.7) 


n  A  a 

from  Eq.  3.6  using  Eq.  3.5. 

For  two  cylinders  in  contact  the  normal  stiffness  at  contacts  is 
based  on  the  Hertz-Mindlin  theory  of  frictional  contacts  with  modifications 
for  conditions  of  local  yielding  at  contacts  (Chang  and  Misra  1988a).  The 


normal  contact  stiffnesses  is  given  by 

1  1  -  v 

K  = 


n  G 


[  2  In  (2R/a)  -  1  j 


(3.8) 


.1/2 


where  a=  (2  (1  -  v)  N  R  /  rt  G)  ,  G  is  an  equivalent  shear  modulus  less 
than  the  elastic  shear  modulus,  \>  is  the  Poisson's  ratio,  R  is  the 
equivalent  radius  of  the  cylinders. 


3.1.2  Non-Hertzian  theory 

The  assumptions  and  restrictions  made  in  the  Hertz  theory  of  elastic 
contact  are  parabolic  profiles,  frictionless  surfaces  and  elastic  half  space 
theory.  However  the  real  situation  is  much  more  complicated.  For  example, 
when  two  non-conforming  elastic  bodies  having  continuous  profiles  are 
pressed  into  contact,  the  pressure  distribution  between  them  is  not 
determined  uniquely  by  the  profiles  of  the  bodies  within  the  contact  area. 
Two  further  conditions  have  to  be  satisfied:  (1)  that  the  interface  should 
not  carry  any  tension  and  (2)  that  the  surfaces  should  not  interfere  outside 
the  contact  area. 

Friction  at  the  interface  of  two  non-conforming  bodies  brought  into 
normal  contact  plays  a  part  only  if  the  elastic  constants  of  the  two 
materials  are  different.  The  mutual  contact  pressure  produces  tangential 


I 

» 

I 

j 

displacements  at  the  interface  as  well  as  normal  compression.  If  the  j 

materials  of  the  two  solids  are  dissimilar,  the  tangential  displacements  , 

will,  in  general  be  different  so  that  slip  will  take  place.  Such  slip  will 
be  opposed  by  friction  and  may,  to  some  extent,  be  prevented.  Therefore  a 

t 

central  region  where  the  surfaces  stick  together  and  regions  of  slip  towards 
the  edge  of  contact  may  be  expected.  If  the  coefficient  of  limiting 
friction  was  sufficiently  high  slip  may  be  prevented  entirely. 

When  the  limits  of  elastic  behavior  have  been  exceeded  and  plastic  ^ 

i 

flow  has  begun,  the  plastic  zone  is  fully  contained  by  the  surrounding  \ 

material  which  is  still  elastic.  For  bodies  having  smooth  profiles,  e.g. 
cylinders  or  spheres,  the  plastic  zone  lies  beneath  the  surface.  In  these 
circumstances  the  material  displaced  is  accommodated  by  an  elastic  expansion 
of  the  surrounding  solid.  The  resulting  deformations  are  of  approximately 
the  same  magnitude  as  the  elastic  deformations.  However  the  results  of  the 
Hertzian  theory  needs  to  be  corrected  for  the  yielding. 

As  the  load  is  increased  further,  the  plastic  zone  breaks  out  to  the 
free  surface  and  the  displaced  material  is  free  to  escape  by  plastic  flow  to 
the  sides.  This  is  the  'unconstrained'  mode  of  deformation  analysed  by  the 
theory  of  rigid-plastic  solids.  The  plastic  zone  would  be  expected  to  break 
out  to  the  surface  and  the  unconstrained  mode  become  possible  when  the 
pressure  at  the  contact  reaches  the  value  given  by  the  rigid-plastic  theory. 

The  contact  behavior  acounting  for  these  Hertzian  restrictions  and 
assumptions  have  been  primarily  studied  by  numerical  techniques.  No  simple 
closed-form  solutions  of  contact  deformations  or  stiffness  for  non-Hertzian 
contacts  exit  in  lietrature  to  the  authors  knowledge. 


3.2  CYCLIC  TANGENTIAL  LOADING  ON  CONTACT 


23 


For  a  constant  normal  force,  the  effect  of  a  tangential  force  T  <  uN 
is  to  cause  'slip'  over  part  of  the  contact  area.  Slip  is  initiated  at  the 
circumference  of  the  circular  contact  area  and,  as  T  increases,  an  annular 
area  of  slip  develops  spreading  radially  inwards  until,  when  T  =  uN,  rigid 
body  sliding  occurs.  If  the  tangential  force  is  reduced  then  the  slip  in 
the  opposite  direction  spreads  radially  inwards  from  the  perimeter  of  the 
contact  area.  Thus,  all  load  reversals  cause  slip  reversals  that  propagates 
radially  inwards  from  the  perimeter  of  the  contact  area,  instead  of  receding 
the  existing  slip  annulus.  Consequently  the  tangential  stiffness  is 
dependent  on  the  loading  history. 

The  tangential  stiffness  is  also  dependent  on  the  magnitude  of  the 
normal  force.  Fig.  3.1  shows  schematically,  the  loading,  unloading  and 
reloading  behavior  for  values  of  normal  force  at  the  contact.  Mindlin  and 
Deresiewicz  (1953)  have  identified  several  loading  sequences  involving 
variations  of  both  normal  and  tangential  forces  and  developed  theoretical 
solutions  for  the  tangential  stiffness. 

A  general  procedure  based  on  Mindlin  and  Deresiewicz  (1953)  work  has 
been  reported  by  Thornton  and  Randall  (1987).  The  procedure  is  to  update 
the  normal  force  and  contact  area  radius  followed  by  calculating  AT  using 
the  new  values  of  N  and  a.  By  reanalysing  the  loading  cases  considered  by 
Mindlin  and  Deresiewicz  (1953)  it  can  be  shown  that  for  loading,  unloading 
and  reloading,  the  tangential  incremental  displacement  may  be  expressed  as 


66  =  8G  a  [“ 
except  when 


U  AN  + 


U  AN 
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and  sign  (A6)  =  sign  (AN) 


(3.9) 


(3.10) 


From  Eq.  3.9 


DISPLACEMENT  -  d 


where 


(3.12) 

(loading)  (3.13) 

(unloading)  (3.14) 

(reloading)  (3.15) 

and  the  negative  sign  in  Eq.  3.11  is  only  evoked  during  unloading,  T  and 
T  are  defined  in  Fig  3.1.  For  a  current  state  given  by  point  1  in  Fig  3.1 
(during  loading,  unloading,  or  reloading)  a  tangential  incremental 
displacement  corresponding  to 

|A&|  =  8  vith  m  >  0  (3.16) 

will  result  in  a  new  state  given  by  point  2  on  the  curve  corresponding  to 
the  new  value  of  N.  Larger  values  of  |A6|  will  result  in  a  state  farther 
along  the  curve  such  as  point  3.  A  problem  occurs  if  the  conditions  given 
in  Eq.  3.10  is  satisfied,  since  point  2  is  not  reached  and  the  new  state 
does  not  lie  on  the  curve  corresponding  to  the  new  value  of  N.  This  case 
can  be  solved  by  setting  0=1  in  Eq.  3-11  until  8GaIA5  >  pEAN. 

For  two  cylinders  in  contact  the  shear  stiffness  at  contacts  based  on 
the  Hertz-Mindlin  theory  of  frictional  contacts  with  modifications  for 
conditions  of  local  yielding  at  contacts  (Mindlin  1949,  Chang  and  Misra 
1988a)  is  given  by 

■£"  -  “T  Y  P  ‘  T/(tan  \  N)]‘1/2  (3.17) 

s  n 

where  ♦  is  the  friction  angle  between  the  two  particles,  N  and  T  are  the 
normal  and  shear  force  at  the  contact  respectively,  and  t  is  a  constant. 


(2  -  vx)  (2  -  v2) 


3  1  -  (T  +  u AN) 

M  N 

1  -  (T*  -  T  +  2uAN) 

9  =  — TTn - 

-31-  (T  -  T**+  2uAN) 
2  Vi  N 
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Chapter  4 

FABRIC  TENSOR  AND  STIFFNESS  PROPERTIES  OF  RANDOM  PACKINGS 


Among  the  significant  geometrical  features  that  influence  the  stress- 
strain  behavior  of  a  granular  assembly  are  the  void  ratio,  the  coordination 
number  and  the  spatial  distribution  of  the  vectors  joining  the  centroid  of 
the  particle  to  the  contact  point  on  the  surface  called  the  contact  position 
vectors.  For  convenience  the  distribution  function  of  the  contact  position 
vectors  can  be  represented  in  a  tensorial  form.  The  distribution  function 
is  related  to  the  'fabric  tensor'  of  the  granular  packing.  Using  this 
distribution  function  and  assuming  a  contact  force-displacement  law  the 
stiffness  constants  for  anisotropic  packings  of  spheres  and  discs  can  be 
obtained  (Chang  and  Misra  1988c,  1988d). 

4.1  GEOMETRICAL  PROPERTIES  OF  THE  PACKING  STRUCTURE 

Geometrical  properties  of  granular  packings  have  been  of  interest  in 
various  areas  of  engineering.  Various  efforts  have  been  made  to  identify 
the  geometrical  characteristics  of  importance  in  the  study  of  granular  media 
in  general  (Gray  1968,  Shahinpoor  1983,  Oda  et.  al.  1982).  Studies  suggest 
that  among  the  geometrical  properties  of  granular  packings  that  influence 
the  mechanical  behavior  are  the  average  coordination  number,  the  void  ratio, 
and  the  distribution  of  the  contact  position  vectors  in  space. 

4.1.1  Coordination  number  and  void  ratio 

For  a  granular  assembly,  the  coordination  number  is  defined  as  the 
number  of  contact  points  per  particle.  The  average  coordination  number  of  a 
random  assembly  is  a  useful  measure  of  the  closeness  of  the  packing  or  the 


LTVTWnJ^ 


a 


void  ratio.  The  void  ratio  is  defined  as  the  ratio  of  the  volume  of  voids 
with  the  volume  of  solids  in  a  packing.  The  total  number  of  contacts  N  in  a 
given  volume  V  of  a  granular  packing  can  be  obtained  from  the  knowledge  of 
the  void  ratio  e,  the  coordination  number  in  and  the  particle  size.  For  a 
packing  of  equal  sized  particles,  the  ratio  of  number  of  contacts  to  the 


volume  of  the  packing,  y,  can  be  expressed  as  follows 


N  3m  „ ,  ,  „  N  m 

=  - - -  (for  spheres);  y  =  - 5 - 

4 n  pJ(l  +  e)  it  p  (1  +  e) 

where  p  is  the  particle  radius. 


(for  discs) 


(4.1) 


Further,  the  correlation  between  the  coordination  number  and  the  void 
ratio  of  the  packing  have  been  studied  experimentally  by  several 
investigators  (Oda  et.  al.  1982,  Yanagisawa  1983).  The  relationship  between 
m  and  e  have  been  experimentally  found  for  lead  and  glass  balls  by  Smith, 
Foote  and  Busang  (1929),  Filep  (1936),  Field  (1963),  and  Oda  (1977). 
Experimental  results  are  also  obtained  by  Marsal  (1973)  for  rockfill 
materials  and  Yanagisawa  (1983)  for  gravels  with  round  and  flat  shapes. 


Some  of  the  empirical  equations  are  listed  as  follows: 

1  v - j -  /mim  i  =  3.  i83^2,469-e^ 


1.  Yanagisawa  (1983) 


2.  Smith,  et.al.  (1929) 


3.  Filep  (1936) 


4.  Field(1963) 


in  =  26.486 


10.726 
(1  -  n) 


(1  +  e) 


(4.2) 


(4.3) 


(4.4) 


(4.5) 


These  four  empirical  equations  are  plotted  in  Fig.  4.1  and  compared 
with  experimental  data.  Experimental  results  from  Oda  et.  al.  (1982) 


v.v.v.v 
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suggests  that  the  relation  between  the  coordination  number  and  the  void 
ratio  is  independent  of  grain  size  distribution. 

4.1.2  Fabric  tensor 

Another  important  measure  of  the  packing  structure  of  granular  media 
is  a  second  order  symmetric  tensor,  termed  the  'fabric  tensor',  that  has 
been  introduced  by  some  investigators  (Ref.  Satake  1982,  Oda  et.  al.  1982, 
Cowin  1985).  More  general  form  of  the  'fabric  tensor'  has  been  introduced 
by  Kanatani  (1984).  In  general,  the  'fabric  tensor'  represents  the  average 
of  the  tensor  product  of  the  observed  directional  data,  such  as  contact 
position  vectors.  The  contact  position  vector  is  defined  as  the  vector 
joining  the  centroid  of  a  particle  to  a  contact  point  (representing  contact 
with  a  neighboring  particle)  on  its  surface  (r11  and  rP  in  Fig.  2.2).  For 
packings  of  spherical  particles  the  contact  position  vector  is  same  as  the 
contact  normal  vector.  If  only  the  unit  vector  n  in  the  direction  of  the 
contact  normal  vector  is  considered  then,  in  general,  the  'fabric  tensor' 
can  be  defined  as  a  tensor  of  rank  m  given  by 

N.  .  =  [  n.  n.  ...n.  S(y,6)  dS  =  <n.  n.  ...n  >  (4.6a) 

where  <  >  denotes  the  ensemble  average  of  the  m  observed  contact  normals 


n.  ,  n.  ...,  n.  ,  t(y»(3)  is  the  density  function  of  the  directional 
11  l2  m 

distribution  of  the  contact  normals,  and  da  =  siny  dy  d(3  is  the  differential 
solid  angle  associated  with  a  unit  sphere  with  0<y<n  and  0 <0<2 n  (y  and  8  are 
defined  in  Fig.  2.3),  such  that, 


r2"  .Ji 

N  .  =*  n.  n.  ...n  S(y,6)  siny  dy  d6 

1  2’*  m  J0  J0  1  2  xm 


(4.6b) 


The  'fabric  tensor'  defined  here  is  useful  in  defining  the  tensorial  form  of 
the  distribution  density  function  £(y,6).  The  relationship  between  the 


'fabric  tensor'  and  the  density  function  £(y,0)  will  be  discussed  in  the 
next  section. 


4.1.3  Distribution  density  of  contact  normals 

The  empirical  distribution  of  the  contact  normals  in  space,  for  a 
random  granular  system,  is  in  general  expected  to  be  discrete  and  random. 

For  continuum  modelling  of  granular  packings,  it  is  useful  to  represent  this 
empirical  discrete  distribution  of  contact  normals  by  a  smooth  density 
function.  Such  a  directional  distribution,  in  three  dimensions,  can  be 
expressed  as  a  spherical  harmonics  expansion  given  by 


CO 

*<Y’P)  •  ^  [l  ♦  L'{  ak0  Pk(cosY) 


k  m  \-i 

P^(cosy)  [a^cos  mg  +  bkmsin  mel  (A‘7) 


m=l 

where  y  and  3  are  defined  in  Fig.  2.3.  Here,  L  denotes  summation  with 
respect  to  even  indices  only,  P^(cosy)  is  the  kth  Legendre  polynomial, 
P^(cosy)  is  the  associated  Legendre  function  and  a^Q,  a^,  and  bkm  are 
parameters.  In  order  to  ensure  that  the  density  function  £(y,0)  is  centro- 
symmetric,  i.e.  E(y>0)  =  S(  h+y,  d+ix) ,  only  the  even  harmonics  are  admissible. 
It  is  evident  that  the  first  term,  i.e.  1,  in  the  expansion,  Eq.  4.7, 
represents  a  sphere  and  the  subsequent  terms  can  be  regarded  as  a  function 
defined  on  the  surface  of  the  sphere.  Further,  since  the  Legendre 
polynomials  and  the  associated  Legendre  functions  are  orthogonal  to  1,  it 
follows  that  the  Eq.  4.7  satisfies  the  identity 

[  E(y,  3)  dQ  =  1  (4.8) 

J  n 


signifying  it  is  a  density  function. 

Similarly,  in  two  dimensions,  the  directional  distribution  can  be 
expanded  in  form  of  a  Fourier  series  given  by 
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t(0)  *  jit  1  +  X]  a^cos  k9  +  b^.  sin  k9  j 


(4.9) 


Due  to  the  centro-symmetric  nature  of  the  distribution,  that  is  £,(9)  = 
t(9+rt),  only  the  even  terms  appear.  Clearly,  Eq.  4.9  also  satisfies  the 


identity  [  £(9)  d9  =  1. 

J  a 


The  expression  for  the  distribution  of  contact  normals,  i.e.  Eq.  4.7 


or  4.9,  can  be  alternatively  written  as  a  cartesian  tensor  equation 

£(n)  *  7—  [l  +  4.,n.n.  +  * .  . ,  ,  n .  n  .  n,  n,  +  . 1 

4n  L  ij  i  J  ljkl  i  j  k  1  J 


(4.10a) 


*(n)  =  ^  [l  *  ♦ijn.n.  +  ♦ijklninjnkn1  «•  . ]  (4.10b) 

which  represents  a  polynomial  in  terms  of  direction  cosines  of  n.  In  Eq. 

4.10,  4^  is  a  coefficient  tensor  of  rank  m  of  appropriate  choice  such 

1’  '  '  m 

that  Eq.  4.10  expands  to  Eq.  4.7  or  Eq.  4.9.  It  is  evident  that  the 
coefficient  tensor  ^  .  can  be  expressed  in  terms  of  the  coefficients  of 

Eq.  4.7  or  4.9. 

Since  Eq.  4.7  represents  a  spherical  harmonics  expansion,  Kanatani 

(1984)  has  shown  that  the  coefficient  tensor  4^  ^  is  a  traceless  tensor. 

1 '  '  ’  m 

He  has  further  shown  that  the  coefficient  tensor  4.  .  can  be  determined 

i. ...  i 

1  m 

using  the  'fabric  tensor'  N.  and  hence  is  directly  related  to  the 

1  ■%  •  •  •  1 
1  m 

packing  structure  of  the  assembly.  The  'fabric  tensor'  N.  .  is 

1l12,**1m 

related  to  the  coefficient  tensor  4.  .  .  .  For  example,  for  m=2  the 

1112"  -1m 

coefficient  tensor  4..  is  related  to  'fabric  tensor'  N.  .  as 

*ij  =  T-  [Nij  -  J  *ij  ]  (for  spheres); 

4..  =  4  I"  N..  -  =■  5..  1  (for  discs)  (4.11) 

ij  L  i)  2  13  J 


U  '.I  t  .<>  ‘ 


Mathematical  details  on  the  spherical  harmonics  expansion  are  discussed  in 
(Butkov  1968,  Kanatani  1984). 


4.2  STIFFNESS  CONSTANTS  FOR  RANDOM  PACKINGS 

Using  the  general  expression  of  the  constitutive  tensor  for  packings 
of  given  by  Eq.  2.41,  the  density  function  of  the  form  given  by  Eq.  4.7  or 
4.9  and  assuming  a  local  constitutive  lav,  the  stiffness  constants  for  a 
random  assembly  of  any  packing  structure  can  be  obtained. 

4.2.1  Constant  contact  stiffness 

Closed  form  solutions  for  the  stiffness  constants  are  derived  by 

assuming  that  all  contacts  in  the  assembly  have  the  same  mechanical  property 

independent  of  the  stress  state.  A  simple  diagonal  form  of  the  local 

constitutive  matrix  relating  the  relative  displacement  vector  d ^  , 

d  ,  d„ )  and  the  incremental  force  vector  f.=(f  ,  f  ,  f  )  (or  d.=(d  ,  d  )  and 
s  t'  l  n’  s  t'  3  n’  s 

f.=(f  ,  f  )  in  two  dimensions)  is  assumed,  such  that,  K11=  K  ,  K„„=K--,=K 
ins  ll  n  zz  s 

and  K. .=  0  (i/j).  K  and  K  are  the  normal  and  the  shear  stiffness 
ij  ns 

respectively  at  the  contact.  It  is  assumed  that  there  is  no  coupling 
between  the  shear  and  the  normal  forces  at  the  contact  such  that  the  off 
diagonal  terms  of  the  local  constitutive  matrix  are  zero  (K^=0  for  i*j). 

4 . 2 . 1 . 1  Packings  of  equal  spheres 

For  the  simplicity  of  further  derivation,  a  truncated  form  of  the 
expansion  in  Eq.  4.7.  consisting  of  second  order  terms  is  used.  The 
Legendre  polynomial  of  degree  2,  i.e.  P-,(cosy),  is  given  by  (Abramovitz  and 


Stegun  1965) 

P2(cosy)  =  1  (3cos^y  -1) 


(4.12) 


The  associated  Legendre  function  P^cosy)  can  be  obtained  from  the  Rodrigues 
formula  (Abramovitz  and  Stegun  1965) 


ii'l’tUVi  Jl’i 


For  k=2  and  m=2r  Eq.  4.13  yields 

2  2 
P2(cosy)  =  3  sin  y 


(4.14) 


Thus  the  truncated  expansion  is  given  as 
$(y.P)  =  4^  f1  +  \  a20  (3  cos2y  +  1)  + 

3  sin^y  (a22C0s2g  +  b22sin2g)  j  (4.15) 

Alternatively,  writing  Eq.  4.15  as  a  cartesian  tensor  equation 
too  -  [l  ♦  ♦ijVj  ]  (4.  If 

where  n  =  (cosy,  siny  cosg,  siny  sing),  and  the  coefficient  tensor  is 
given  by 


(4.16) 


[*..]=  0 

ij* 


2  a20  +  3a22 


(4.17) 


2  a20  "  3a22 


The  first  term  in  Eq.  4.16,  clearly,  represents  the  isotropic  portion  of  the 
distribution  while  the  second  term  represents  the  anisotropic  part.  For 
brevity,  the  isotropic  and  the  anistropic  parts  can  be  combined  as  follows 

«”)  ■  4l  [<5ij  *  V  Vj  ]  ■  4n  [Fij  "i"j  ]  <4'18> 

where  F^  is  a  generalized  'fabric  tensor'  that  characterizes  the 

distribution  function  of  the  normal  vectors  at  the  contacts.  Substituting 

Eq.  4.18  into  Eq.  2.41,  the  constitutive  tensor  becomes 


■f  Jo” ",  th  k 


"lqpk  “  4n  V  JO  JO  "q  Ali  "ij  Tjk  np  Frs  Vs  siny  dy  df3 

Taking  into  account  the  symmetry  of  the  incremental  stress  and  strain 

tensors,  the  stress-strain  relationship,  equivalent  to  Eq.  4.19,  can  be 

expressed  in  Voigt's  notation  as  follows, 

<«  >  =  Cmn  <fen>  (4.20) 

m  mn  n  ' 

where 


(4.19) 


<b  >  =  C  <e  > 
m  mn  n 
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(4.32) 


C45=  -WT-  L  —  (Kn  +  fV  J  (A*33) 

The  other  elements  of  the  stiffness  matrix  C  are  zeros.  The  ratio  of 

mn 

N 

number  of  contacts  to  the  volume  of  the  packing,  — ,  is  expressed  in  terms  of 
the  coordination  number  m  and  void  ratio  e  by  Eq.  4.1. 

4. 2. 1.2  Packings  of  equal  disks 

For  simplicity,  a  truncated  form  of  the  Fourier  expansion  in  Eq.  4.9 
vith  terms  upto  fourth  order  are  used.  The  truncated  expansion  is  given  as 


(4.35) 


£(9)  =  Jn  [3  +  a2  cos^®  +  ^2  +  a4  cos49  +  sin49  j  (4.34) 

Alternatively,  Eq.  4.34  can  be  written  as  a  cartesian  tensor  equation 
«(n)  -  ji  [l  .  ♦ijVj  ♦  ♦,.kln.n.Vl  ]  (4.35) 

where  n  =  (cos©,  sin9). 

Taking  into  account  the  symmetry  of  the  incremental  stress  and  strain 
tensors,  the  stress-strain  relationship,  equivalent  to  Eq.  4.19,  can  be 
expressed  in  Voigt's  notation  as  follows, 


<b  >  =  C  <e  > 
m  mn  n 

where  <b  >  =  (b  ,  b  ,  b  )  and  <s  >  =  (e  ,  e  ,  2c  ). 

m  yy  xx  yx  n  yy  xx  yx 

The  stiffness  constants  C  are  given  as 

mn  6 

cii=  a.n,r7)  (6K  +  2K  )  +  4a,  K  +  a.  (K  -  K  ) 

11  8n  (1+e)  L  n  s  2n  4  n  s  J 


(4.36) 


C,,=  - 7  (6K  +  2K  )  -  4a-  K  + 

22  8n  (1+e)  L  n  s'  2  n 


a.  (K  -  K  ) 
4  n  s 


(4.37) 


(4.38) 
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(4.39) 


'33“  16n  (1+e) 


4  (K  +  K  )  -  2a,  (K  -  K  ) 
n  s  4  n  s 


(4.40) 


C, 


13“  16 n  (1+e) 


2b-  K  +  b.  (K  -  K  ) 
2  n  4  n  s 


(4.41) 


C-.,=  - r  2b-  (K  +  K  )  -  b.  K 

23  16n  (1+e)  L  2  '  n  s'  4  n 

It  is  noted  that  for  the  case  of  isotropic  packing  structure,  i.e. 
a2=a^=b2=b^=0,  Eqs.  4.37-4.42  yield  similar  results  as  Bathurst  and 
Rothenburg  (1987). 


(4.42) 


4.2.2  Hertzian-Mindlin  Contact 

Under  the  condition  for  which  the  shear  force  does  not  exceed  the 

frictional  strength  at  a  contact,  the  local  constitutive  law  can  also  be 

defined  based  on  frictional  Hertzian-Mindlin  contact.  The  contact 

stiffnesses,  K  and  K  ,  obtained  from  the  Hertz-Mindlin  theory  for 
ns 

frictional  contacts  (Mindlin  and  Deresiewicz  1953)  are  given  in  Chapter  3. 

If,  instead  of  a  linear  local  constitutive  law,  a  non-linear 

Hertzian-Mindlin  contact  law  is  assumed,  the  expressions  of  C  become  very 

mn 

complex  and  the  close  form  solutions  can  only  be  obtained  for  the  case  of 
statistically  isotropic  packing  under  initial  isotropic  stress  condition. 

For  a  statistically  isotropic  packing,  the  contact  force  vector  f, 
under  initial  isotropic  stress  (<7^=  is  the  same  for  a11  contacts, 

given  by 


Thus  using  Eqs.  4.21-4.33  with  a2Q=a22='322=^’  t*ie  ^°^sson,s  ratio  and  the 


Young's  modulus  for  the  packing  are  given  as 
v 


v. 


iso  “  2  (5-3v) 
2 


_  4  r  (5-4v) 

^iso  30  (5-3v) 
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f 

For  other  packings  and  initial  loading  conditions  other 


(4.44) 

(4.45) 

than  isotropic,  the 


expressions  become  very  complex  and  a  numerical  approach  to  the  solutions  is 
more  tractable. 


4.3  TYPES  OF  MATERIAL  SYMMETRY 

Due  to  the  different  packing  structure,  a  granular  assembly  may, 

conceivably,  exhibit  different  types  of  material  symmetries.  The  forms  of 

the  stiffness  matrix  for  various  types  of  material  symmetries  are  shown  in 

Table  4.1  (Nye  1957).  Evidently,  the  various  types  of  material  symmetries 

can  be  represented  by  an  appropriate  choice  of  the  'fabric  tensor'  F_. 

Along  with  the  forms  of  stiffness  matrix  for  the  different  material 

symmetries  the  corresponding  'fabric  tensor'  F„  are  presented  in  Table  1. 

The  packing  structure  of  a  granular  packing  with  orthotropic  material 

symmetry  non-coincident  with  the  stress  axes  can  be  represented  by  a  second 

rank  symmetric  fabric  tensor  with  three  non-zero  parameters  a^Q,  a22>  and 

t>22‘  A  granular  packing  in  this  case  will  have  13  independent  non-zero 

stiffness  constants  given  by  Eqs.  4.21  to  4.33.  It  is  noted  that  for  this 

material  the  normal  stresses  a  ,  a  ,  <r  are  coupled  with  the  shear  strain 

z  z  yy  xx 

e  . 
xy 

Granular  packings  with  orthotropic  material  symmetry  can  be 
represented  by  a  diagonal  second  rank  fabric  tensor  with  two  non-zero 
parameters  &2q  an^  a22'  Orthotropic  material  is  characterized  by  symmetry 


TABLE  -  1  For as  of  stiffness  matrices  C  and  the  fabric  tensor  P. . 

an  ij 

for  various  material  symmetries 


Stiffness  Matrix  C 


Fabric  tensor  F 


Ortho tropic  vit 
non-coincident  stress  axes 
‘  *  *  *■ 

*  *  .  * 

*  * 

symm  *  * 

* 

* 


1  +  a„ 


2  a20  +  3a22 


2  a20  ~  3a22 


Tetragonal 

r  *  * - * 

* 


1  +  a20 

0 

o 

0 

1  +  a20 

0 

0 

0 

1  -  2a2( 

Transverse  isotropic 
r  *  * - *  . 

* 


1  +  a„ 


2  a20 


1  2  a20 
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about  three  mutually  orthogonal  planes.  The  stiffness  matrix  of  materials 
with  orthotropic  symmetry  is  represented  by  9  independent  constants.  For 
orthotropic  materials  the  normal  stresses  are  not  coupled  with  the  shear 
strains  and  the  vice  versa. 

Granular  packings  with  tetragonal  material  symmetry  have  6 
independent  stiffness  constants.  Tetragonal  symmetry  is  characterized  by  a 
an  axis  of  symmetry.  If  the  parameters  defining  the  fabric  tensor  F„  were 
chosen  such  that  a2Q  =  2 a22  and  b22=0,  then  packings  with  tetragonal 
symmetry  can  be  represented. 

In  addition  to  possessing  an  axis  of  symmetry,  certain  granular 
assemblies,  have  an  additional  condition  on  the  stiffness  constants,  that  is 
«  j  ^22  “  C23^'  num^er  independent  constants  characterizing  the 

stiffness  matrix  for  such  packings  are  thus  reduced  to  5.  This  type  of 
material  symmetry  is  termed  transverse  isotropic  symmetry.  The  fabric 
tensor  for  the  transverse  isotropic  material  can  be  represented  by  choosing 

a22=b22=0' 

If  the  fabric  tensor  is  chosen  such  that  a20=a22=b22=0’  then  packings 
with  isotropic  material  properties  can  be  represented.  It  is  noted  that  the 
stiffness  constants  for  packings  with  different  material  symmetry  can  be 
obtained  from  Eqs.  4.21  to  4.33  by  using  appropriate  values  of  a2Q,  a22  and 


22' 


4 . 4  DISCUSSION 

To  demonstrate  the  effect  of  packing  anisotropy  on  the  mechanical 
properties  of  random  packings,  parametric  studies  are  performed  for  the 
statistically  isotropic,  transverse  isotropic  and  orthotropic  packings. 

Fig.  4.2  shows  the  density  functions  £(n)  of  the  contact  normal  vectors  for 


y  -  0  plant 


x  ■  0  plant 


z  ■  0  plant 


Fig.  4.2  Section  of  density  function  E(n)  for  statistically  transverse 
isotropic,  orthotropic,  and  monoclinic  packings  along  x=0,  y=0, 
and  2=0  planes. 


the  monoclinic,  the  orthotropic  and  the  transverse  isotropic  packings. 
Three  sections  of  the  density  function,  along  x=0,  y=0,  and  2=0  planes 
respectively,  are  plotted  in  Fig.  4.2.  The  following  values  of  the 


parameters  are  used  for  this  plot:  3~2q~  0.8,  a22=^2Z=^  ^or  transverse 

isotropic  packing,  &2q=  0.2,  a22=0-25,  ^22=0  for  orthotropic  packing,  and 

&2o=  0.2,  322=0.25,  b22=0.1  for  monoclinic  packing. 

The  solutions  are  computed  based  on  Eqs.  4.21  to  4.33  using  the 

following  values  for  the  parameters:  n  =  9,  e  =  0.57,  r=0.01  in.,  and  Kn= 

3800  psi.  Fig.  4.3  shows  the  Young's  modulus  and  the  Poisson's  ratio 

plotted  against  the  ratio  K  /K  (the  shear  stiffness  K  to  the  normal 
r  °  s  n  s 

stiffness  K^)  for  statistically  isotropic  packing.  As  expected,  the  Young's 

modulus  increases  and  the  Poisson's  ratio  decreases  with  the  ratio  K  /K  . 

s  n 

The  maximum  value  for  the  Poisson's  ratio  is  0.25  when  the  ratio  K  /K  is 

s  n 

zero  for  perfectly  smooth  spheres  with  no  tangential  interactions. 

Figs.  4.4,  4.5  and  4.6  show  the  variation  of  the  Young's  moduli,  the 
Poisson's  ratios,  and  the  shear  moduli,  respectively,  with  the  ratio  Ks/Kn 
for  the  transverse  isotropic  and  the  orthotropic  packings.  The  Young's 


moduli  E  =E  and  E  ,  the  two  Poisson's  ratios  v  =v  and  v  ,  and  the 
xx  yy  zz’  zx  zy  xy 

shear  moduli  G  =G  and  G  for  the  transverse  isotropic  are  computed  for 
zx  zy  xy 

the  parameter  a2Q=0.8.  It  is  noted  that  the  shear  modulus  G^^  can  be 
expressed  in  terms  of  Exx  and  v  ,  thus  there  are  only  5  independent 
constants  for  the  transverse  isotropic  packing.  The  three  Young's  moduli 


E  ,  E  and  E  ,  the  three  Poisson's  ratio  v  ,  v  and  v  ,  and  the  three 
xx  yy  zz  zx  zy  xy’ 

shear  moduli  G  ,  G  and  G  for  the  orthotropic  packings  are  plotted  for 
zx  zy  xy 

the  value  of  the  parameter  a2Q=0.2  and  parameter  a22=0-25. 

To  show  the  effect  of  degree  of  anisotropy  on  t he  mechanical 
properties,  the  ratios  of  Young's  moduli  for  transverse  isotropic  packing  to 


RATIO 


Fig.  A. 5  Variation  of  Poisson's  ra 
statistically  transverse 
respectively. 


SHEAR  MODULUS  (psixlO5) 


RATIO  K,/Kn 


0.0  0.2  0.4  0.6  0.8  1.0 

RATIO  Ka/Kn 

4.6  Variation  of  Shear  moduli  with  the  ratio  for  statistically 

transverse  isotropic  and  orthotropic  packings  respectively- 


the  Young's  modulus  for  isotropic  packing  are  plotted  against  the  parameter 


in  Fig.  4.7.  The  ratios  of  the  Poisson's  ratios  for  the  two  packings 
also  shown  in  Fig.  4.7.  Clearly  the  more  anisotropic  packing  structure 
a  more  anisotropic  mechanical  properties. 

5 

I 

t 

£ 

» 

I 

s 

U 


a20 

are 

has 


V 

V 
'■i. 


Chapter  5 

EXPERIMENTAL  STUDY  ON  ROD  ASSEMBLIES 

Experimental  tests  vere  performed  to  provide  a  database  for  the 
verification  of  the  developed  theory.  Experiments  vere  performed  to  obtain 
stress-strain  response  of  packings  with  regular  and  random  structures. 
Specifically  of  interest  was  the  directional  dependence  of  the  stiffness 
properties  of  the  packings.  Details  of  the  experimental  setup  are  presented 
in  Xue  1988. 

5 . 1  APPARATUS 

The  apparatus  used  in  this  study  consists  of  a  loading  frame  composed 
of  4  rigid  aluminium  plates  1/2"  thick  and  2"  high.  The  maximum  size  of  rods 
assembly  that  can  be  placed  in  the  loading  frame  is  8"  by  8"  in  dimension. 
Each  rod  of  assembly  is  placed  vertically,  resting  on  a  glass  table  top 
(shown  in  Fig.  5.1).  The  frame  is  so  designed  that  it  can  be  compressed 
independently  in  X  and  Y  directions.  The  box  can  also  be  made  to  swing  about 
Y-axis,  hence  it  can  serve  as  a  directional  shear  box. 

The  loading  is  applied  using  a  set  of  push-pull  type  pneumatic 
cylinders  (or  pistons)  controlled  by  regulators  with  air  pressure  gages 
attached.  Cylinders  can  be  pulled  or  pushed  to  apply  compressive  forces  and 
shear  forces  to  this  box  with  desired  stress  conditions,  such  as:  isotropic 
compression,  biaxial  compression,  and  shear  loading  with  a  controlled 
rotation  of  principal  stress  axis.  Five  dial  gages  that  read  to  0.0001"  are 
attached  (see  Fig.  5.1)  to  monitor  deformations. 


5.2  SPECIMENS  AND  MEASUREMENTS 


Six  regular  and  two  random  packings  have  been  investigated.  Samples 
of  approximately  300  aluminum  cylinders,  of  2  in.  high,  are  placed  randomly 
or  in  accordance  with  a  given  regular  packing  arrangement. 

Regular  packing  structres  selected  in  this  study  consist  of  four 
packings  of  equal  sized  cylinders  and  two  packings  of  multi-sized  cylinders 
as  shown  in  Fig.  5.2.  For  packing  A,  B,  C  and  D,  the  diameter  of  the 
cylinders  used  is  0.5  inch.  For  packing  E,  the  diameters  of  the  cylinders 
used  are  0.5  inch  and  0.25  inch,  and  for  packing  F,  0.75  inch  and  0.5  inch. 
Two  types  of  random  packings  are  used,  namely  one  sized  packing  made  of  271 
rods  of  0.5  in.  diameter  and  two  sized  packing  made  of  229  rods  of  0.5  in. 
diameter  and  190  rods  of  0.25  in.  diameter. 

Strains  of  the  rods  assembly  is  obtained  from  the  deformation  of  the 
box  measured  from  the  attached  5  dial  gages.  The  thickness  of  the  frame 
plates  is  designed  to  be  large  enough  such  that  the  bending  deflection  of 
the  loading  frame  is  insignificant. 

Stresses  applied  on  the  rods  sample  are  calculated  from  the  readings 
of  the  pressure  gages  connected  to  the  pneumatic  cylinders.  Since  the 
actual  forces  transmitted  to  the  rods  assembly  are  different  from  those 
directly  calculated  from  pressure  gages  due  to  friction  loses  in  the  system, 
calibration  was  made  to  account  for  the  friction  loses  in  the  piston  of  the 
pneumatic  cylinders  and  the  friction  loses  between  the  loading  frame  and  the 
glass  plate.  In  addition,  the  stress  madnitude  applied  to  the  assembly  was 
selected  to  be  much  larger  than  the  friction  between  rods  and  the  smoothened 
glass  plate  so  that  the  effect  of  friction  is  expected  to  be  insignificant. 

Photographs  of  the  assembly  packing  were  taken  for  the  random 
assembly.  The  rods  are  painted  black  on  the  top  and  their  centers  are 


Packing  A 


Packing  C 


Packing  E 


marked  so  that,  from  photographs,  the  coordinates  of  the  center  point  of 
each  circle  can  be  digitizd  to  give  the  contact  normal  distribution. 

5.3  STRESS  PATHS 

The  stress  path  applied  to  the  sample  is  a  shear  test  with  fixed 
principal  stress  axis  as  illustrated  in  Fig.  5.3  on  a  normalized  deviatoric 
stress  vs  shear  stress  plot.  All  tests  start  with  isotropic  initial  stress 
state,  represented  by  the  point  of  origin  in  Fig.  5.3.  An  initial  isotropic 
stress  of  2.5  psi  is  used  for  all  the  tests.  Then  the  incremental  stress  of 
0.5  psi  is  applied  to  the  sample  in  a  direction  with  an  angle  a  inclined 
from  Y-axis  such  that  the  principal  stress  axis  is  rotated  by  an  angle  of  a, 
with  the  major  principal  stress  a^=  3.0  psi  and  the  minor  principal  stress 
o^=2. 5  psi.  For  each  packing,  tests  were  conducted  for  the  following  a 

angles  (0°,  15°,  22.5°,  30  0°,  37.5°,  and  45.0°). 

5 . 4  TEST  DEVIATION 

Each  test  was  repeated  several  times  to  evaluate  the  test 

variability.  The  coefficients  of  variations  of  the  measured  strains  t  ,  e 

x  y 

and  shear  strain  e  are  found  to  be  about  20%,  18X  and  15%  for  e  ,  e  and 
xy  x  y 

e  respectively.  The  scatter  may  be  caused  by  many  factors.  In  addition 
xy 

to  the  factors  associated  with  loading  system  and  apparatus,  one  factor 
identified  is  that  the  rods  are  not  perfectly  circular  and  equal  in 
diameter.  Thus  identical  samples  are  not  possible  to  be  reconstructed. 


followed  for  the  tests  on  rod  assembl 


Chapter  6 

COMPARISON  OP  THEORETICAL  AND  EXPERIMENTAL  RESULTS 


Verification  of  the  theory  is  carried  out  by  comparing  the 
theoretical  computations  with  the  measured  experimental  results. 

Comparisons  are  made  with  the  experimental  results  from  the  tests  on  rod 
assemblies  in  directional  shear  box.  Comparisons  are  also  made  with  test 
results  from  glass  balls  assembly  in  a  resonant  column  device  and  tests  on 
sands  under  low  amplitude  waves.  Response  under  different  initial  stress 
conditions  are  also  computed  and  compared  with  measured  trends. 

6.1  COMPARISON  WITH  EXPERIMENTS  ON  ROD  ASSEMBLIES 

6.1.1  Regular  packings 

To  compare  with  the  experimental  results,  the  response  of  the  six 

regular  packings  (shown  Fig.  5.2)  under  the  loading  conditions  discussed 

above  were  computed  using  the  theory.  The  values  of  the  constants  used  for 

2 

theoretical  computations  in  this  work  are:  G  =  15000  psi  (105000  kN/m  ),  f  = 
2.5,  v  =  0.1,  and  ♦  =  15°. 

7  Tp 

The  apparent  deformation  modulus  of  a  packing  is  defined  by 
Aa 

E  =  “i-2  <6-X> 

a 

where  Aa  is  the  incremental  stress  in  the  direction  a,  e  is  the  measured 
a  a 

incremental  normal  strain  in  the  direction  a,  calculated  from  measured 
deformation  of  test  box.  The  apparent  Poisson's  ratio  can  be  given  by 


where  €  is  the  normal  strain  measured  in  the  direction  perpendicular  to 
that  of  £a-  Normal  strains  £a  and  £a+gQ  are  not  necessary  the  principal 
strains. 

Figs.  6.1  and  6.2  show  the  apparent  moduli  and  apparent  Poisson's 
ratio  for  the  six  packings  respectively.  The  measured  results  are  compared 
with  the  predicted  theoretical  results  for  various  values  of  a  on  polar 
coordinates . 

Packing  A  (square  packing)  is  most  stiff  when  loaded  in  the  direction 
of  Y-axis.  When  it  is  loaded  in  the  diagonal  direction  (i.e.,  as  a  becomes 
to  45  degree)  the  modulus  is  minimum  while  the  Poisson's  ratio  is  maximum. 
Packing  B  (hexagonal  packing)  behaves  isotropically.  Both  moduli  and 
Poisson's  ratio  are  equal  in  all  directions.  Packing  C,  is  a  mixture  of 
square  packing  and  hexagonal  packing.  Remarkably  the  modulus  and  Poisson's 
ratio  of  packing  C  are  average  of  those  of  packing  A  and  packing  B.  Packing 
D  shows  the  similar  isotropic  behavior  as  that  of  packing  B,  however,  the 
magnitude  of  the  modulus  of  packing  D  is  much  smaller  than  packing  B.  Both 
packing  E  and  F  are  two-sized  packings.  Their  moduli  and  Poisson's  ratios 
are  anisotropic.  More  detailed  comparisons  are  included  in  Chang  and  Misra 
1988a,  Chang,  Misra  and  Xue  1988. 

6.1.2  Random  packings 

Two  types  of  random  packings  are  used  for  the  experiments.  Random 
packing  A  is  made  of  one  size  rods  with  total  number  of  contacts  N=1372,  the 
void  ratio  e=0.179,  and  the  coordination  number  m=5.06.  Random  packing  B 
consists  of  rods  of  two  different  sizes  with  total  number  of  contact  N=1991, 
the  void  ratio  e=0.179,  and  the  coordination  number  m=4.75.  The 
distribution  of  the  branch  vectors,  vector  joining  the  centroids  of  two 
particles  in  contact,  is  shown  in  Figs.  6.3  and  6.4  for  the  two  packings. 
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Fig.  6.1 


Apparent  moduli  vith  direction  for  the  six  regular  packings. 
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Fig.  6.2  Apparent  Poisson's  ratio  with  direction  for  the  six  regular 
packings . 
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ig.  6.4  Branch  vector  distribution  for  tvo-size  random  packing  of 
cylinders. 
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Clearly  the  branch  vector  is  two  times  the  contact  position  vector  for  the 
case  of  one  size  packings. 

To  compare  with  the  experimental  results,  three  sets  of  values  of  the 
constants  to  represent  stiff,  soft  and  average  response  were  used.  The 
values  of  the  constants  used  are:  G  =  40000  psi  (280000  kN/m^),  f  =  1.0  for 
stiff  packing;  G  =  15000  psi  (105000  kN/m^),  T  =  2.5  for  soft  packing;  and  G 
=  25000  psi  (175000  kN/m^),  ?  =  2.0  for  the  best  fit  of  the  experimental 
results.  For  all  the  computations  the  value  of  \>  is  0.1,  and  is  15°. 

The  theoretical  predictions  and  the  measured  strains  are  shown  in 
Figs.  6.5  and  6.6  plotted  against  the  direction  of  the  principal  stress  axis 
a.  The  measured  data  is  completely  bounded  within  the  predicted  response. 

The  theory  gives  reasonable  predictions  of  the  trend.  Measured  strain 
shows  more  scatter  since  it  is  very  small  in  magnitude  and  therefore 
expected  to  have  more  measurement  errors. 


6.2  COMPARISON  WITH  EXPERIMENTS  ON  GLASS  BALLS  IN  RESONANT  COLUMN  DEVICE 

To  evaluate  the  applicability  of  the  theory,  the  theoretical  results 
are  also  compared  with  the  experimentally  observed  shear  modulus  of  packings 
of  uniform  glass  balls  from  resonant  column  method.  A  confining  pressure 
is  initially  applied  to  the  cylindrical  samples  made  of  glass  balls 
contained  within  a  rubber  membrane  (fanagisawa  1983).  Then  the  torsional 
oscillation  is  applied  on  the  top  of  the  specimen  to  develop  a  shear  wave 
transmitted  through  the  specimen.  The  shear  wave  velocity  or  the  shear 
modulus  of  the  material  can  be  obtained  from  the  observed  resonant  frequency 
and  the  dimensions  of  the  specimen. 

Since  no  information  is  available  on  the  packing  structure,  for 
computing  the  theoretical  results  we  assume  that  the  packing  of  glass  balls 
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Fig.  6.5  Strains  versus  direction  of  major  principal  stress  axis  a  for  the 
one-size  random  packing  of  cylinders. 
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is  statistically  isotropic.  The  shear  modulus  is  taken  to  be  5.75x10^  psi 

7  2 

(4.0x10  kN/m  )  and  the  Poisson's  Ratio  is  0.13  for  glass.  The  relationship 
between  the  co-ordination  number  m  and  the  void  ratio  e  used  in  this 
prediction  is  given  by  Eq.  4.2  by  Yanagisawa  (1983). 

Fig.  6.7  shows  the  comparison  of  the  predicted  and  measured  moduli 
for  the  packing  of  glass  balls.  As  can  be  seen  from  Fig.  6.7,  the  shear 
modulus  is  influenced  significantly  by  the  confining  pressure  a  .  The 
experimental  values  show  that  the  shear  moduli  increase  in  proportion  to  the 
0.42  power  of  the  confining  pressure  <tc  while  in  theory  the  power  is  0.33. 
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Eq.  6.4  is  shown  in  Fig.  6.8  compared  with  measured  data  on  void  ratio  and 
coordination  number.  For  convenience,  a  linear  line  is  used  to  fit  the 
measured  data.  The  comparison  between  the  shear  moduli  obtained  from  the 
theory  and  the  empirical  equation  is  shown  in  Fig.  6.9. 

6.3.2  Effect  of  initial  stress  conditions 

The  initial  stress  condition  is  an  important  factor  influencing  the 
mechanical  behavior  of  a  packing.  A  closed-form  solution  is  difficult  to 
obtain  for  initial  stress  conditions  other  than  isotropic  for  the  reason 
that  the  forces  vary  at  contacts  and  hence  the  contact  stiffness  is 
different  for  each  contact. 

Moduli  of  an  isotropic  packing  are  computed  for  different  initial 
stress  conditions.  The  different  initial  stress  conditions  are  represented 
by  points  of  stress  state  lying  on  three  stress  paths  shown  in  Fig.  6.10 
namely,  1)  compression  path,  2)  constant  mean  path,  and  3)  reduced 

2 

compression  path.  Initially,  stresses  a2=  10  psi  (69  kN/m  ).  In 

path  1,  the  axial  stress  a ^  is  increased  while  the  confining  stress  (o^o^) 
is  kept  constant.  In  path  2,  is  increased  and  a 2  and  (^decreased  while 
the  mean  stress  is  kept  constant.  In  path  3,  and  CT2  are  decreased  while 
the  vertical  stress  is  kept  constant.  The  moduli  are  computed  using 
G=7.5xl0^psi  (5.2x10^  kN/m^),  v=0.13,  and  <^=24°.  It  is  observed  that, 
under  these  stress  conditions,  the  packing  exhibits  cross-isotropic  behavior 
with  properties  symmetrical  about  z  axis. 

The  computed  Young's  moduli  in  the  vertical  direction  are  shown  in 
Fig.  6.11  for  the  three  stress  paths.  With  an  increase  of  stress  ratio,  the 
moduli  in  the  vertical  direction  increase  for  the  compression  path  while  the 
moduli  decrease  for  the  reduced  compression  path.  For  the  constant  mean 
path,  the  moduli  decreases  only  slightly.  This  trend  is  in  agreement  with 
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the  trend  observed  in  the  experimental  tests  by  Yanagisava  (1983).  For 
stress  ratio  greater  than  2.3,  the  Young's  moduli  and  the  shear  moduli 

shov  appreciable  decrease  because  the  shear  force  exceeds  the  contact  shear 
strength  for  many  contacts.  For  comparison,  the  computed  Young's  moduli  in 
the  horizontal  direction  and  the  shear  moduli  are  also  shown  in  Fig.  6.11 
and  6.12  respectively. 
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Fig.  6.12  Shear  moduli,  G  and  G 


versus  stress  ratio  a. /a. 
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